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A key observable in strongly interacting resonant Fermi gases is the contact parameter C, which 
governs both the pair correlation function at short distances and the momentum distribution at 
large momenta. The temperature dependence of C was recently measured at unitarity, where ex- 
isting theoretical predictions differ substantially. We report accurate data for the contact and the 
momentum distribution in the normal phase of the unitary gas, obtained by Bold Diagrammatic 
Monte Carlo. In our scheme, C is extracted from the pair correlation function, while the C/k 4 tail 
of the momentum distribution, being built in at the analytical level, is free of /c-dependent noise. 

PACS numbers: 05.30.Fk, 67.85.Lm, 74.20.Fg 



The resonant Fermi gas is a fundamental model of 
quantum many-body physics. It features a smooth 
crossover between fermionic and bosonic superfluidity, 
as predicted in the context of condensed matter physics 
[1-4 and confirmed by remarkable experiments on ul- 
tracold atomic Fermi gases near Feshbach resonances [5] . 
It is also relevant to neutron matter [6] and high-energy 
physics [7], particularly in the unitary regime reached at 
the center of the crossover. As a result of the vanishing 
interaction range, resonant Fermi gases feature charac- 
teristic ultraviolet singularities governed by an observable 
called contact [4|l8HTT]. In particular, the density-density 
correlation function at short distance diverges as 

and the momentum distribution has the tail 

^ W (2) 

Here C is the contact per unit volume for the homo- 
geneous gas, n a (r) = ^(r)^ a (r) is the density oper- 
ator, and the spin-a momentum distribution n a (k) is 
normalised to j 7v(k)d 3 fc/(27r) 3 = n a = (n a (r)). A di- 
rect manifestation of Eq. ([!]) is that in a unit volume, 
the number of pairs of fermions separated by a distance 
smaller than s is Cs/(A-k) in the s — » limit. Hence 
C controls the (anomalously high) density of pairs with 
vanishing interparticle distance [8] [TUJ [12] . 

The contact is directly related to a large variety of 
observables that have been measured by various exper- 
imental techniques: the population of the closed chan- 
nel molecular state measured by laser molecular spec- 
troscopy p~3j [14] , the large- momentum tail of the static 
structure factor measured by Bragg spectroscopy [T5l - 
[17], the tail of the momentum distribution measured by 



non-interacting time-of-flight or by momentum-resolved 
radiofrequency spectroscopy [18] . the derivative of the 
energy with respect to the inverse scattering length [9 
extracted from the pressure equation of state measured 
by in-situ imaging [19], and the large- frequency tail in 
radiofrequency spectroscopy [18, 20] [2Tj . 

The recent experimental study [2T] , being both lo- 
cally resolved and temperature dependent, is the first to 
have access to the temperature dependence of the con- 
tact of the homogeneous unitary gas. On the theoretical 
side, however, current predictions based on summation 
of certain classes of Feynman diagrams [22H24] or lattice 
Quantum Monte Carlo simulations [25] differ substan- 
tially even at the qualitative level, especially on approach 
to the superfluid transition from the normal side. 

In this Letter, we present controlled accurate results 
for the contact and the momentum distribution of the 
unitary Fermi gas in the normal phase, obtained by Bold 
Diagrammatic Monte Carlo (BDMC). We evaluate the 
contact C by expressing it as a density of pairs. We fur- 
ther demonstrate that the momentum distribution ob- 
tained within our BDMC scheme features the C/k A tail 
with the same value of C and without /c-dependent sta- 
tistical noise. 

In our BDMC scheme, all Feynman diagrams are 
stochastically sampled to high order, and abelian resum- 
mation is applied to the diagrammatic series [26] [27] . The 
diagrams are built on single-particle and pair propagators 
which are bold, i.e. fully dressed. We directly extract the 
contact from the bold pair propagator T thanks to the 
relation 

C = -r(r = 0,r = 0~) (3) 

(we set h and m to unity). This relation is consis- 
tent with the interpretation of C in terms of a density 
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FIG. 1: Resummation of the dimensionless contact CA 4 (at 
fiji — 1). After linear extrapolation to e — >• 0, consistent 
results are obtained by the different resummation methods: 
Lindelof (o), shifted Lindelof (x), and Gauss (+). Also shown 
are the first order [23] and third order results without resum- 
mation (solid and dashed horizontal lines). 



of pairs, since it is formally analogous to the relation 
n a = G a (r = 0, r = 0~) between the single-particle den- 
sity and the bold single-particle propagator G. Although 
Eq. (§ was first obtained within T-matrix approxima- 
tions [24] [28j [29], it is actually exact. A simple way to 
derive it is to use the regularized version of Eq. which 
holds in a lattice model [30] , 



C = g 2 (n t (0) %(0)) 



(4) 



i.e. the contact is equal to the double occupancy, up to a 
renormalization factor set by the bare coupling constant 
go (see also Q2]). The result Q then follows from the 
fact that 



T(r,T)=g S(T) 



b 3 



or diagrammatically 



V (T(^ t )(r,r)(^)(0,0)) 

(5) 



(6) 




where T is the time-ordering product, and the first term 
does not contribute in the continuum limit where the 
lattice spacing b tends to zero. 

We start by verifying the consistency between different 
resummation methods. As in Ref. [26] we use abelian 
resummation methods, where the N-th order terms of 
the series are multiplied by e _eAiV - 1 , with A n depend- 
ing on the resummation method: A n = nlnn for Lin- 
delof, A n = n 2 for Gauss, and A n = (n — 1) ln(n — 1) for 
shifted-Lindelof. Figure [T] shows the e-dependence of the 
dimensionless contact CA 4 , for a fixed value of the grand- 
canonical parameter P/a. Here j3 = l/(/c#T) is the inverse 



temperature and A = \^2tt/3 is the thermal wavelength. 
When the regularization parameter e is extrapolated to 
zero, the three methods give consistent results. To under- 
stand the overall shape of the resummation curves and 
to make a link with the self-consistent T-matrix approx- 
imation of Refs. [23, 31 , we remark that for e — >• oo, the 
Lindelof and Gauss curves would asymptote to the result 
obtained by truncating the diagrammatic series at or- 
der 1 as in Refs. [23, 31 , while the shifted-Lindelof data 
would tend to the substantially different result obtained 
by truncating the diagrammatic series at order 3. 
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FIG. 2: BDMC results for the contact in the high-temperature 
regime, compared to the virial expansion [24] [32] at order two 
(dashed line) and three (solid line). 

The temperature dependence of the contact in the 
high-temperature regime is shown in Fig. [2j where we 
compare with the virial expansion 



CA 4 = 16tt 2 (c 2 



c 3 e 



(?) 



The coefficients c 2 = 1/tt [32] and C3 ~ —0.141 [24] come 
from the 2-body and 3-body problem respectively. The 
agreement with the third-order coefficient C3 constitutes 
a non-trivial test of our approach, since this coefficient 
comes from the three-body propagator [33] which has 
to be reconstructed in our case by summing an infinite 
number of diagrams. 

The behavior of the contact in the low-temperature 
region of the normal phase is displayed in Fig. [3j The 
contact in terms of canonical variables, C(n, T) or equiv- 
alently C/k F 4 versus T/Tp, is shown in Fig. [3^1. Our 
results are consistent within error bars with the experi- 
mental results of 21] . The lattice Auxiliary Field Quan- 
tum Monte Carlo data of Ref. [25 disagree with our data 
and contain a priori unknown systematic errors com- 
ing from discretization of space (i.e. finite filling factor) 
and of imaginary time. In contrast with the non self- 
consistent T-matrix results of Ref. [22 and the Nozieres 
Schmitt-Rink results of Ref. [24] , we do not observe any 
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FIG. 3: Temperature dependence of the contact in the low-temperature region of the normal phase, in terms of (a) canonical 
and (b) grand canonical variables. BDMC (this work): blue solid circles, JILA experiment [21]: brown solid diamonds, lattice 
AFQMC simulations [25] : grey crosses. The curves correspond to different diagrammatic approximations: non self-consistent 
T-matrix [22] : dashed black line, self-consistent T-matrix [23] : solid red line, Nozieres-Schmitt-Rink [24]: dotted green line in 
(a) and dashed black line in (b). The transition point to the superfluid phase from Refs. [34, 35] is indicated by the green 
arrow. 



pronounced enhancement when decreasing temperature. 
The self-consistent T-matrix results of Ref. [23] are re- 
markably close to our data in Fig. [3^. The contact in 
terms of grand canonical variables, C(/i,T) or equiva- 
lently C//u 2 versus fin, is shown in Fig. [3)3. It is nat- 
ural to use these variables to discuss the different dia- 
grammatic results since the diagrammatic technique is 
formulated in the grand-canonical ensemble. The func- 
tions C(n, T) and C(/i,T) are connected by the equation 
of state n(/i,T), given for each of the considered ap- 
proaches in Refs. [26, 36-38 . The non self-consistent T- 
matrix and Nozieres-Schmitt-Rink approaches yield the 
same result for CY/i,T). Indeed, in both approaches, C 
is given by Eq. (I3| with T replaced by T^ ) which is the 
sum of the ladder diagrams built on the ideal gas propa- 
gators . The self-consistent T-matrix approach over- 
estimates C(/i,T) as seen in Fig. j3|), and also overesti- 
mates the equation of state n(/i,T) (see [26]). These two 
systematic errors largely cancel out when one considers 
C(n,T) as in Fig.^. 

We now turn to the momentum distribution. The 
large-momentum physics governed by the contact en- 
ters in the following way the bold-line diagrammatic 
formalism (see also [24, 39). In short, the tail of 
the momentum distribution comes from the diagram 
depicted in Fig. [4j More precisely, for k — » 00 the 
Dyson equation simplifies to n a (k) = G^k, 0~) ~ 

1% dn f% dr 2 G ( °\k, -n)E a (k, n - T 2 )G ( °\k, r 2 ). 

Here the large-momentum ideal gas propagators G^ can 
be replaced by vacuum propagators, which are retarded 




FIG. 4: Leading diagrammatic contribution to the momen- 
tum distribution n C7 (k) at large /c, which can be interpreted 
physically as the simultaneous propagation of two opposite- 
spin particles of large and nearly opposite momenta and of a 
missing pair of lower momentum p <C k. Imaginary time runs 
from right to left. The single-particle lines propagate forward 
in time and can be replaced by vacuum propagators. The pair 
propagator runs backwards in time and is fully dressed. 



and quickly decaying functions of imaginary time. Hence 
the integral is dominated by T2 — » + and t\ — > 0~, and 
the imaginary time argument t\ — of the self-energy 
tends to 0~. The behavior of S(k, r) for k — >• 00, r — » 0~ 
can be determined analytically by considering the bold 
diagrammatic series. The lowest-order diagram reads 



/^r(p,T)G_ ff (-k+p,-T), where 



S (1) (k,r) 

the integral is dominated by momenta p <C k and the 
propagator G can be replaced by the vacuum expression 
_ e -efc|r| i n relevant region 0<— t<1/& 2 — ^0 
(here e& = k 2 /2). As for the higher order bold diagrams, 
their contribution is suppressed because they include 
integrations over internal times which are restricted to 
narrow ranges, since G and T become narrow functions 
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of time at large momentum. Using Eq. (|3| we thus 
obtain 

E(k, r) ~ C e efcT , k oo, r -> 0". (8) 

After substituting this form into the asymptotic Dyson 
equation given above, the large- momentum tail, Eq. ([2|, 
is recovered. 

These analytical considerations have the following im- 
plications for our BDMC calculation. Since the C/k 4 tail 
of the momentum distribution comes exclusively from the 
lowest-order self-energy diagram, this tail is automati- 
cally built into our self-consistent BDMC scheme pro- 
vided this diagram is evaluated with high precision. We 
achieve this by using numerical Fourier transformations 
(rather than Monte Carlo) and analytical treatments of 
leading-order singularities [27 , in the spirit of [31 . Our 
value of C, of course, still differs from the one of the self- 
consistent T-matrix approximation of Refs. [23j[3l], since 
our lowest-order diagram is built on the fully dressed pair 
propagator T given by the BDMC self-consistency which 
includes higher-orders contributions. On the technical 
side, we note that treating the lowest-order self-energy 
diagram separately (without using Monte Carlo) has an- 
other advantage: the steep function of r in Eq. (|8| would 
be hard to capture by Monte Carlo sampling. 

The momentum distribution multiplied by k 4 is shown 
in Fig. [5] for our lowest temperature. The large- 
momentum tail is reproduced without /c-dependent sta- 
tistical noise (in contrast to other Monte Carlo meth- 
ods [25j HQ]) and perfectly agrees with our value of the 
contact determined from Eq. (|3|. The momentum distri- 
bution for various temperatures is shown in Fig. [6| As 
the temperature is lowered, we do not observe any in- 
cipient discontinuity at fop, as would be the case for a 
pronounced degenerate Fermi liquid behavior. 

In conclusion, we have incorporated the ultraviolet 
physics governed by the contact into the framework of 
BDMC. We accurately evaluated the contact via the bold 
pair propagator. Moreover, our momentum distribution 
data automatically contain the large-momentum tail with 
the proper contact coefficient, as we have shown analyt- 
ically and confirmed numerically. Our data can serve as 
benchmarks for future experimental and theoretical stud- 
ies, not only on the contact and momentum distribution 
themselves, but also in various contexts where the con- 
tact appears either in ultraviolet asymptotics [23j I4TH44] 
or in sum rules [HI [42j [45j [46] . We expect that the possi- 
bility to incorporate analytical knowledge into BDMC is 
a generic capability of the approach which will facilitate 
its applications to other outstanding strongly correlated 
fermion systems. 
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FIG. 5: BDMC data for the momentum distribution n a (k), 
multiplied by k 4 in order to reveal the large-momentum tail 
n a (k) — »> C/k 4 . Dashed line: value of the contact C com- 
puted directly from the pair propagator. The temperature is 
T/T F = 0.191(5) (/3/i = 2.25). 
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FIG. 6: BDMC data for the momentum distribution at vari- 
ous temperatures. 
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